Correlated sampling in quantum Monte Carlo: a route to forces 
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In order to find the equilibrium geometries of molecules and solids and to perform ab initio 
molecular dynamics, it is necessary to calculate the forces on the nuclei. We present a correlated 
sampling method to efficiently calculate numerical forces and potential energy surfaces in diffusion 
Monte Carlo. It employs a novel coordinate transformation, earlier used in variational Monte Carlo, 
to greatly reduce the statistical error. Results are presented for first-row diatomic molecules. 
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Over the past decade, quantum Monte Carlo (QMC) 
methods have been used to calculate the struc- 

tural and electronic properties of a variety of atoms, 
clusters and solids. For systems with large numbers of 
electrons, QMC methods at present provide the most 
accurate benchmark calculations of structural energies. 
However, a major difficulty of QMC methods has been 
the determination of equilibrium geometries and poten- 
tial energy surfaces. Hence, most QMC calculations have 
been performed on geometries obtained with either den- 
sity functional theory (DFT) or conventional quantum 
chemistry methods. The computation of forces on nu- 
clei has been a stumbling block that has limited a more 
widespread use of QMC methods. 

DFT methods or standard quantum chemistry tech- 
niques use the Hellman-Feynman theorem to compute 
the forces on nuclei [Q. Unfortunately, this is not prac- 
tical within QMC for three reasons. First, the wave 
functions used in QMC calculations are usually not ob- 
tained by minimizing the energy. Therefore, if Hellman- 
Feynman theorem were employed in variational Monte 
Carlo (VMC), the forces would have a systematic error. 
Second, in fixed-node diffusion Monte Carlo (DMC) the 
Hellman-Feynman force has an error due to the disconti- 
nuity in the derivative of the fixed-node wave function at 
nodes Q. Finally, in both VMC and DMC, the statisti- 
cal fluctuations would be too large since the fluctuations 
of the potential energy are much larger than those of the 
total energy. 

Alternatively, one could simply compute energy differ- 
ences to obtain either forces (for an infinitesimal displace- 
ment of the ions) or the full potential energy surface of 
the system. However, while quantum chemistry meth- 
ods can rely on having an approximately constant and 
smoothly varying error in the energy, a major disadvan- 
tage of QMC methods is that, in addition to systematic 
errors, one has statistical errors which make the determi- 
nation of energy differences or smooth potential energy 
surfaces very expensive in computer time. Even though it 



is not possible to entirely eliminate the statistical errors, 
it is possible, by using correlated sampling j^, to make 
the statistical errors in the relative energies of different 
geometries much smaller than the errors in the separate 
energies and to make them vanish in the limit that the 
two geometries become identical. In the past, the corre- 
lated sampling technique has been used within VMC 
but there have been very few attempts joj to extend the 
approach to DMC, and these were approximate and/or 
inefficient and were tested only on H2, Hj^ and LiH. 

In this paper, we present a novel DMC correlated sam- 
pling technique to efficiently compute accurate forces and 
potential energy surfaces. The DMC bond lengths of 
first-row diatomic molecules computed with this algo- 
rithm are found to be in better agreement with experi- 
ment values than are the VMC, Hartree-Fock (HF), local 
density approximation (LDA) and generalised gradient 
approximation (GGA) values. 

Correlated sampling in variational Monte Carlo. In- 
stead of performing independent VMC runs which would 
have independent statistical errors, one generates MC 
configurations for a reference situation only. The MC 
configurations are sampled from where ip is the wave 
function for the reference situation. Then, unbiased ex- 
pectation values for somewhat different secondary wave 
functions V's are obtained by reweighting the configura- 
tions sampled from ■0^, e.g., for Hamiltonians Ti. and Tig, 
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where the weights of the A'conf MC configurations are 



W,. = 
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and R = (ri, . . . , r^f). The statistical error in E^ — E is 
considerably smaller than that in E^ or E, making this 
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method of calculating Eg — E more efficient than per- 
forming independent VMC computations of Es and E. 

Space-warp coordinate transformation. Since we want 
to compute the relative energy of two different geome- 
tries, the Hamiltonians Tl and Hs in Eq. |^ correspond 
to two sets of nuclear coordinates, and R^, for the 
reference and the secondary geometry, respectively. 

The electronic coordinates sampled from the reference 
wave function ip'^ will not be optimal for computing the 
energy Eg corresponding to the nuclear coordinates R^, 
since the electron density will be peaked at Rq, rather 
than at R^. To solve this problem, a mapping of the 
electron coordinates was introduced in Ref. Q such that 
those electrons that are close to a given nucleus move 
almost rigidly with that nucleus: 
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(We use Latin indices for electronic coordinates and 
Greek indices for nuclear coordinates.) F{r) is any suffi- 
ciently rapidly decaying function such as r^^, exp(— Kr) 
or exp(K/r). The reduction in statistical error due to 
the use of these warped electronic coordinates for the sec- 
ondary geometries is large and almost independent of the 
choice for F{r) and k. Unless otherwise stated, all results 
in this paper were computed with F{r) = r~'^ and k — A. 
The equation for E^ — E (Eq. P is now 
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Ef=r|V's(Rp/^(R,)|'^(R. 
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and J(R) is the Jacobian for the transformation (Eq. |^). 

Correlated sampling in diffusion Monte Carlo. In 
DMC, the imaginary-time evolution operator exp(— Tir) 
is used to project out the ground state from the trial 
wave function within the fixed-node and the short-time 
approximations [lO| ]. The primary walk is generated ac- 
cording to the stochastic implementation of the integral 
equation: 



/(R',i + T) = ydRG(R',R,r)/(R,i) 



(7) 



where G'(R',R,r) = (R'| exp{-HT}|R). If we introduce 
importance sampling using the reference wave function ij), 



the importance sampled Green's function for small values 
of r (short-time approximation) is given by the product 
of three factors, diffusion, drift and growth/decay: 

G(R',R,r) = -—^e-' ^^e^(^ ■^■^), (8) 

(zttt) 2 

where V = V?/'(R)/?A(R) and S'(R',R,t) = {E^ - 
£;l(R') - -Bl(R))t/2 with E^L = HV'(R)/V'(R-)- At time 
t, a set of primary walkers characterized by the pairs 
{R,i^Wi) is a random realization of the distribution /: 



/(R,t) = ^w, J(R-R0. 



(9) 



Each walker executes a branching random walk: a walker 
originally at R drifts to R + V(R)t and then diffuses to 
R' according to the Gaussian term in Eq. ||. To ensure 
that, in the limit of perfect importance sampling (i.e. V 
is the ground state wave function ip^), we are sampling 
ijp' despite the short-time approximation in the Green's 
function, the move is accepted with probability 



p = min < 1, 



|V>(ROpr(R,R',r) 
|^(R)|2f(R',R,T) 
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as prescribed by the detailed balance condition. We de- 
note by T the drift-diffusion part of the Green's func- 
tion G. Finally, the weight of the walker is multiplied 
by exp[S'(R', R, r)]. In practice, we employ the modi- 
fied version of G presented in Ref. |jll|] that takes into 
account the non-analyticities of V and i?L at the nodes 
and particle coalescence points. 

Given a primary walk generated according to Eq. ^, the 
secondary walk is specified by the space-warp transfor- 
mation (Eq. H). Two complications, absent in VMC, arise 
for correlated sampling in DMC. First of all, the dynam- 
ics of the secondary walker should have been governed 
by an importance sampled Green's function constructed 
from the secondary wave function -03, G's(R''', R*', r), and 
the move should have been accepted with probability 



Ps = min { 1 



|7^,(R-0pf,(R%R-',T) 
|7^,(R'^)|2f;(R«',R^r) 
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However, the secondary-geometry move was effectively 
proposed according to the drift-diffusion Green's func- 
tion T(R',R, rVJ(R') and accepted with probabilty p 
defined in Eq. To correct for the wrong dynamics, we 
should multiply the weights of the secondary walkers by 



a(R-^R^T) 

T(R',R, r)/J(R') 



r — 



(12) 



where r — Ps/p if the move is accepted and r = (1 — 
Ps)/(1 ~ p) if the move is rejected. However, these prod- 
ucts fluctuate wildly (r can be anywhere between zero 
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and infinity). Therefore, it is not practical to follow this 
route to perform correlated sampling unless bounds can 
be placed on the ratios while at the same time ensuring 
that unbiased results are obtained in the r — > limit. 

An additional complication is the common practice in 
fixed-node DMC to reject moves that cross nodes. If a 
primary walker attempts to cross, the move is rejected, 
p in Eq. |l^ is set to zero, and the ratio r for the sec- 
ondary walker (Eq. ^2|) becomes ill-defined. Moreover, if 
primary and secondary walkers were to be treated on the 
same footing {ps set to zero when the secondary walker 
crosses its own nodes), the weights of the secondary walk- 
ers would all become zero in a sufficiently long run. Even 
though this problem can be easily overcome since it is le- 
gitimate to do fixed-node DMC allowing walkers to cross 
nodes |pl| , reweighting as in Eq. ^ remains impractical 
due to the large fluctuations. 

In this paper, we propose an alternative correlated 
sampling algorithm. Our algorithm is approximate but 
very accurate. Given the successful implementation of 
correlated sampling within VMC, we wish to devise a 
scheme that differs as little as possible from VMC but 
yields results very close to the fixed-node DMC limit 
for the secondary geometries. We perform the pri- 
mary walk as described above and generate the sec- 
ondary walks according to the space-warp transforma- 
tion. In the averages, we retain the ratios of the sec- 
ondary and primary wave functions as would be done in 
VMC (Eqs. I and |). The secondary weights are the 
primary ones multiplied by the product of the factors 
exp[5s(R"',R",r^) - 5(R',R,r)] for the last N^.^i gen- 
erations. Note that, in the exponential factors, we intro- 
duced Ts, a time-step for the secondary path, in general 
different from the time-step r used for the primary walk. 
Due to the warp transformation, the secondary moves 
were effectively proposed with a different time-step, Ts, 
in the drift-diffusion term. A sensible definition of r= is 
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where Ai? is the displacement resulting from diffusion, 
and Ai?s is the displacement needed to take the sec- 
ondary walker from its drifted position to the position 
specified by the space- warp transformation. Having com- 
puted Ts over the first equilibration blocks of our DMC 
run, we will use this time-step for computing the drift 
and reweighting of the secondary walk. Therefore, each 
secondary geometry has a different time-step Tg. 

Secondary geometry wave functions. We considered 
three choices for secondary geometry wave functions: 
(1) The secondary geometry wave functions have the 
same parameters {p} as the primary one but the co- 
ordinates are relative to the new nuclear positions: 
V's(Ri,Ra) = V'(Ri,Ra,Ps) with Ps = p, possibly with 
the minimal changes required to impose the cusp condi- 
tions. 



(2) The secondary geometry wave functions at warped 
electron positions are related to the primary ones at the 
original positions, -0s(R^,R^) = '(/;(Ri, R^, p)/-v/7(R~). 
This wave function depends on the transformation (it 
was used in Ref. h-c\ with a different transformation) 
and has the advantage that the weights Wi in (Eq. |^) 
are unity. Surprisingly, it gives larger fluctuations of the 
energy differences than choice (1). 

(3) i/'s(Rj, Rq) — '0(Rj7 Rq7 Ps) with reoptimized param- 
eters Ps. This choice gives the smallest fluctuation of the 
energy differences and the best potential energy surface. 
We calculate all molecules with choice (1) but also 
demonstrate the superiority of choice (3) for B2. 

Results and conclusions. The algorithms presented in 
the previous sections are tested on first-row homonuclear 
dimers. The primary wave functions were optimized 
close to the experimental bond length by the variance 
minimization method The potential energy curves 
were obtained with correlated sampling from ten geome- 
tries, using the warp transformation and recentered sec- 
ondary geometry wavefunctions (choice (1) above). 
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FIG. 1. VMC fluctuations ((Jvmc) of the relative energy of 
the primary and secondary geometries divided by the bond 
stretch for B2 . The smallest avMC is achieved by using warp- 
ing along with reoptimized secondary wave functions. 

To ascertain the efficiency of our method, we per- 
formed two additional calculations for B2; in the first, 
we omitted the warp transformation, whereas in the sec- 
ond we employed reoptimized, rather than recentered, 
secondary wavefunctions. In Fig. [l], we present the VMC 
root-mean-square fiuctuations (uvmc) of the relative en- 
ergy of primary and secondary geometries divided by the 
atomic displacement, AE/AR, for B2. Introducing the 
warp transformation yields a reduction of about a factor 
of 3.5-5 in ctvmCi which corresponds to a factor of 12-25 
saving in computer time. Moreover, (Tvmc is only slightly 
dependent on the secondary geometry used. As expected, 
a further reduction in uvmc is obtained when the space- 
warp transformation is used in combination with reopti- 
mized, rather than recentered, secondary geometry wave 
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functions. The space-warp transformation was found to 
be of even greater help for heavier molecules, e.g. for F2 
the reduction in the fluctuations was a factor of 3.2-8.2. 

To test the accuracy of our DMC correlated sampling 
algorithm, we performed DMC runs for H2 and B2 for 
three different primary geometries, (a) the equilibrium 
geometry, (b) a geometry stretched by 0.2 and (c) by 
—0.2 a.u. The runs (a), (b) and (c) should give identical 
potential energy curves if the algorithm were exact. In 
Fig. H, we show results for B2 that reveal the high ac- 
curacy of our DMC algorithm: the three DMC curves 
are very close and clearly distinguishable from the VMC 
results. These results are confirmed by the calculations 
for H2 where, despite the use of an intentionally poor 
wave function, the three curves gave the equilibrium bond 
lengths (a) 1.4014(2) (b) 1.4014(2) and (c) 1.4015(2) a.u. 
The true equilibrium bond length, from a careful fit to 
the results of Rcf. ||l2|, is 1.4011 a.u. 
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FIG. 2. Potential energy curve for B2 in VMC and DMC. 
The three DMC curves axe obtained with three different pri- 
mary geometries (equilibrium, stretched by 0.2 and —0.2 a.u.) 
and using recentered wave functions. All curves are shifted 
with the energy at the equilibrium distance (arrow) defined 
as the zero. Atomic units are used. 

To test the improvement resulting from employing 
Ts^T (Eq. |l^), we performed, for II2, DMC correlated 
sampling with Tg = r. Since Tg > r for Ai? > and Tg < r 
for Ai? <0, we expect this potential energy curve to yield 
an equilibrium bond-length that is too short. The equi- 
librium bond length is indeed 1.4003(2) a.u., which is 4 
standard deviations from the true bond-length, whereas 
that obtained with our r^^r algorithm, 1.4014(2) a.u., 
is 1.5 standard deviation from the true bond- length. 

Having ascertained the accuracy and efhciency of our 
algorithm, we computed the bond lengths of all first-row 
dimers with VMC and DMC correlated sampling. In Ta- 
ble ^, we list the errors in the bond lengths obtained from 
restricted Hartree-Fock (RHF) LDA GGA 
VMC and DMC. The RHF results show the worst agree- 
ment with experiment, with Be2 not being bound. The 



DMC errors are, in all cases, either smaller than or com- 
parable to those from VMC, and are smaller than LDA 
and GGA errors by a factor of 3.9 and 2.6, respectively. 

In this letter, we presented an efficient method to com- 
pute numerical forces in DMC, a long-standing unsolved 
problem in QMC techniques. The method is very accu- 
rate and was tested on first-row dimers, where the DMC 
bond lengths were found to agree with experiment better 
than those from HF, LDA, GGA and VMC. 



TABLE I. Experimental @ bond lengths (in a.u.) of 
first-row dimers and theoretical errors in RHF, LDA, GGA, 
VMC and DMC. 



molecule 


Expt. 


RHF 


LDA 


GGA 


VMC 


DMC 


Li2 


5.051 


0.270 


0.069 


0.057 


0.101(2) 


0.018(3) 


Bea 


4.630 




-0.109 


-0.001 


-0.069(3) 


-0.014(5) 


B2 


3.005 


0.086 


0.025 


0.042 


0.018(2) 


0.002(2) 


C2 


2.348 


-0.007 


0.006 


0.023 


0.006(2) 


0.008(1) 


N2 


2.074 


-0.061 


-0.006 


0.011 


0.012(2) 


0.007(1) 


O2 


2.282 


-0.107 


-0.012 


0.044 


0.028(2) 


0.023(4) 


F2 


2.668 


-0.161 


-0.053 


0.040 


0.021(4) 


0.015(5) 


rms 




00 


0.054 


0.036 


0.049 


0.014 
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